Experimental enforcement of sex-limited autosome inheritance does not reveal intralocus sexual conflict

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = FALSE)

Load packages

library(Matrix) # my pc needs this to load the tidyverse correctly
library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(loo) # to use information criteria in brms models
library(dagitty) # to build DAGs
library(ggdag) # to tidy DAGs
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables

Supplementary methods

\(~\)

Table S1. Recipe for food medium used in our experiment. The provided quantities make ~ 1 litre of food.

tibble("Ingredients" = c("Soy flour", "Cornmeal", "Yeast", "Dextrose", "Agar", "Water", "Tegosept", "Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)"),
       "Quantity" = c("20 g", "73 g", "35 g", "75 g", "6 g", "1000 mL", "17 mL", "14 mL")) %>% 
  pander(split.cell = 40, split.table = Inf)

Ingredients	Quantity
Soy flour	20 g
Cornmeal	73 g
Yeast	35 g
Dextrose	75 g
Agar	6 g
Water	1000 mL
Tegosept	17 mL
Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)	14 mL

\(~\)

Load in the data

\(~\)

fitness_data <- read_csv("data/SLC_fitness_data.csv") %>% 
  mutate(Fitness_vial_ID = as.factor(Fitness_vial_ID),
         Block = as.factor(Block),
         Population = as.factor(Population),
         Treatment = as.factor(Treatment),
         GFP = as.factor(GFP),
         Sex = as.factor(Sex),
         Rearing_vial = as.factor(Rearing_vial),
         Total_red_offspring = Red_female_offspring + Red_male_offspring,
         Total_bw_offspring = bw_female_offspring + bw_male_offspring,
         Total_offspring = Total_red_offspring + Total_bw_offspring) %>% 
  rename(Inheritance_treatment = Treatment)

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'fitness_data')),
      pageLength = 692
    )
  )
}


my_data_table(fitness_data %>% select(-Comment))

Column explanations

Fitness_vial_ID: a unique identifier for each trial of the fitness assay.

Block: the experiment was run in three distinct blocks, using flies from separate generations.

Population: we measured the fitness of flies from 12 independent populations that contained autosomes that had undergone experimental evolution.

Inheritance_treatment: the populations carried autosomes that had been exposed to one of three inheritance treatments for 20 generations: a female-limited inheritance treatment where autosomes were always passed from mother to daughter, a male-limited treatment where autosomes were passed from father to son, and a control condition where inheritance was unconstrained.

GFP: the GFP marker carried by the population. UBI indicates the presence of a transgene that encodes ubiquitous expression of GFP, while 3xP indicates the presence of a different transgene that encodes the expression of GFP in the ocelli.

Sex: the sex of the individuals that we were measuring the fitness of.

Rearing_vial: the vial the treatment flies used in the trial developed in. This variable is included to capture variation explained by the rearing environment e.g. small differences in food moisture content or quantity. Note that females and males can have the same rearing vial as the sexes were reared together.

Red_female_offspring: the number of adult female offspring sired/produced by flies sourced from one of the 12 populations.

Red_male_offspring: the number of adult male offspring sired/produced by flies sourced from one of the 12 populations.

bw_female_offspring: the number of adult female offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

bw_male_offspring: the number of adult male offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

Total_red_offspring: the total number (sexes pooled) of adult offspring sired/produced by flies sourced from one of the 12 populations.

Total_bw_offspring: the total number (sexes pooled) of adult offspring sired/produced by competitor flies.

Total_offspring: the total number (sexes and eye colours pooled) of adult offspring counted in each vial.

\(~\)

Modelling approach

\(~\)

Female and male fitness are fundamentally different concepts / traits. There are also several differences between our female and male fitness assays. The major difference is that the male assay contains half the number of females in any given vial than does the female assay. The logic behind this design choice is that sexually selected processes are a more important determinant of male fitness than they are female fitness, so any fitness differences may only be observed when competition for fertilisations is high.

For these reasons, we choose to split the data up and model female and male fitness separately.

female_fitness <- 
  fitness_data %>%
  filter(Sex == "Female")

male_fitness <- 
  fitness_data %>%
  filter(Sex == "Male") %>% 
  mutate(prop_red = Total_red_offspring / Total_offspring)

We fit the following fixed and random effects to model female and male fitness. Our aim is to estimate the causal effect that Inheritance_treatment (I) has on fitness (F).

The DAG below shows our understanding of the scientific question we are attempting to model. Arrows indicate the direction of causal effects, at least as we interpret them.

# Lets make a function to speed up the DAG making process

gg_simple_dag <- function(d) {
  
  d %>% 
    ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
    geom_dag_point(color = met.brewer("Hiroshige")[4]) +
    geom_dag_text(color = met.brewer("Hiroshige")[7]) +
    geom_dag_edges() + 
    theme_dag()
  
}

dag_coords <-
  tibble(name = c("I", "G", "B", "P", "R", "F"),
         x    = c(1, 1.5, 2.5, 1.5, 2.5, 2),
         y    = c(2, 3, 3, 1, 1, 2))

dagify(F ~ I + G + B + P + R,
       coords = dag_coords) %>%
  gg_simple_dag()

Fixed effects

Inheritance_treatment (I): this is the inheritance regime that the autosomes carried by each of the populations were subject to for 20 generations. There are three levels: populations carrying female-adapted autosomes, populations containing male-adapted autosomes and populations carrying control autosomes that experienced an unmanipulated inheritance regime. We are designing our model to test for a causal effect of this variable.

Mediator variables

Block (B): fitness might differ between the three distinct blocks we split our experiment up into. Blocks differed temporally, used flies from different generations and different batches of food. It is also possible that there were minor fluctuations in the lighting and temperature environment experienced during development between blocks. Each of these variables may introduce variation into our fitness measurements, that can be accounted for by including the Block variable in our model.

GFP (G): it is possible that fitness may be affected by the GFP transgene carried by each population. For example, one could imagine that any unintended fitness effects of a transgene might be of greater magnitude if it is expressed in a larger proportion of tissues, as is the case for the UBI transgene versus the 3xP transgene. Note that each GFP type is carried by an equal number of populations from each of the three evolutionary treatments.

Varying/Random effects

Population (P): our design contained 12 independent populations of autosomes that originated from a single outbred laboratory fly population. The populations were split and autosomes from each were subjected to one of the three evolution treatments for 20 generations. 4 populations experienced a female-limited inheritance regime, 4 a male-limited regime and 4 an unlimited or control regime.

Rearing_vial (R): the vial individual flies developed within may introduce further variation into our response variable. Like Block this variable controls for micro-environmental variation.

\(~\)

Accounting for over-dispersion

The data is over-dispersed with several highly influential (outlier) observations that have large effects on our posterior prediction. To combat this, we fit models following the betabinomial distribution family, as this is better equipped to deal with extreme values at the tails i.e. overdispersion.

However, the beta-binomial is not a native family in brms, we need to create the distribution family using the custom_family() function. The code below is taken directly from the brms_customfamilies vignette, which can be viewed here.

beta_binomial2 <- custom_family(
  "beta_binomial2", dpars = c("mu", "phi"),
  links = c("logit", "log"), lb = c(NA, 2), 
  # note that we set the lower bound to 2, following McElreath, rather than Buerkner. This means that the most conservative estimate for phi we get is a flat expectation between 0 and 1
  type = "int", vars = "vint1[n]"
)


stan_funs <- "
  real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
    return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
  }
  int beta_binomial2_rng(real mu, real phi, int T) {
    return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
  }
"

stanvars <- stanvar(scode = stan_funs, block = "functions")

\(~\)

Female fitness

\(~\)

To estimate the fitness of females carrying each of the three autosome types, we placed three experimental females into a yeasted vial with three female competitors that carried the bw mutation. We then introduced six males that also carried the bw mutation. We allowed them to mate and oviposit for three days, then removed all adults from the vial. 12 days later we counted all of the adult progeny in the vial and scored them for eye-colour. Progeny with red eyes were produced by the experimental females ( bw is recessive) and progeny with brown eyes were produced by the competitor females. We calculated fitness as the proportion of red eyed offspring in the vial.

We present the fixed effects from the model:

\(~\)

female_fitness_model <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Inheritance_treatment + Block + GFP + (1|Population) + (1|Rearing_vial),
      data = female_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = phi)),
      iter = 4000, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.95, max_treedepth = 10),
      seed = 2, stanvars = stanvars, file = "Fits/female_fitness.model")

fixef(female_fitness_model) %>% 
  kable(digits = 3) %>% 
  kable_styling()

	Estimate	Est.Error	Q2.5	Q97.5
Intercept	1.112	0.095	0.920	1.300
Inheritance_treatmentFemale_limited	0.168	0.104	-0.050	0.376
Inheritance_treatmentMale_limited	0.185	0.103	-0.018	0.388
Block2	-0.637	0.081	-0.793	-0.477
Block3	-0.624	0.072	-0.765	-0.483
GFPUBI	-0.322	0.085	-0.489	-0.154

We need to write some additional code to get some post processing stuff i.e. LOO to work. Code courtesy of the brms_customfamilies vignette, which can be viewed here.

Run LOO to see if we’ve effectively modelled over dispersion.

female_fitness_model <- add_criterion(female_fitness_model, criterion = "loo", file = "Fits/female_fitness.model")

loo(female_fitness_model)

## 
## Computed from 8000 by 327 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1399.9 10.6
## p_loo        23.0  2.1
## looic      2799.8 21.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     326   99.7%   2845      
##  (0.5, 0.7]   (ok)         1    0.3%   1130      
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

The beta-binomial model looks good. It returns no points with high pareto k values.

Conduct a posterior predictive check to confirm our model is doing what we want it to.

pp_check(female_fitness_model, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

The posterior predictive distribution matches our raw data quite well, indicating the model is functioning as we wanted.

\(~\)

Derive predictions from the posterior

Get predictions from the model and present them in a Table

draws <- as_draws_df(female_fitness_model)

draws_female <-
  draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(`Female-limited`, Control, `Male-limited`) %>% 
    pivot_longer(cols = c(`Female-limited`, Control, `Male-limited`),
                 names_to = "Inheritance_treatment") %>% 
  rename(prop_focal_offspring = value) %>% 
  arrange(Inheritance_treatment)

draws_diff <- draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>%
  mutate(`Female - Control` = `Female-limited` - Control,
         `Male - Control` = `Male-limited` - Control,
         `Female - Male` = `Female-limited` - `Male-limited`) %>% 
  select(`Female - Control`, `Male - Control`, `Female - Male`)

 #draws_diff %>% 
  # summarise(`Estimated prop. of offspring produced` = mean(`Male - Control`)*100,
   #         `2.5%` = quantile(`Male - Control`, probs = 0.025) *100,
    #        `97.5%` = quantile(`Male - Control`, probs = 0.975) *100)

draws_female %>% 
  group_by(Inheritance_treatment) %>% 
  summarise(`Estimated prop. of offspring produced` = mean(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance_treatment	Estimated prop. of offspring produced	2.5%	97.5%
Control	0.752	0.715	0.786
Female-limited	0.782	0.747	0.812
Male-limited	0.785	0.753	0.814

# now find the differences between the control and the sex-limited treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

prop_table_female <-
  draws %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Inheritance_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Inheritance_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring produced relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = `difference comparison`) #%>% 
#pander(split.cell = 40, split.table = Inf, round = 3)

\(~\)

Male fitness

\(~\)

To estimate the fitness of males carrying each of the three chromosome types, we conducted an experiment very similar to the female fitness assay. However, because male fitness has stronger covariance with fertilisation events than does female fitness, we conducted the male fitness assay with a 1:2 sex ratio (female:male) rather than the 1:1 ratio used in the female assay. This increases the strength of sexual selection and is a more appropriate way to expose differences in fitness between groups of males. As with the females, we calculated fitness as the proportion of red-eyed offspring in the vial.

male_fitness_model <-
  brm(Total_red_offspring | vint(Total_offspring) ~ 1 + Inheritance_treatment + Block + GFP + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = beta_binomial2, 
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = phi)),
      iter = 4000, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.95, max_treedepth = 10),
      seed = 2, stanvars = stanvars, file = "Fits/male_fitness.model")


fixef(male_fitness_model) %>% 
  kable(digits = 3) %>% 
  kable_styling()

	Estimate	Est.Error	Q2.5	Q97.5
Intercept	0.765	0.233	0.299	1.221
Inheritance_treatmentFemale_limited	0.417	0.268	-0.133	0.946
Inheritance_treatmentMale_limited	0.105	0.266	-0.430	0.631
Block2	0.970	0.149	0.670	1.265
Block3	-0.034	0.142	-0.304	0.247
GFPUBI	-0.336	0.226	-0.786	0.117

Run LOO to see if we’ve effectively modelled over dispersion.

male_fitness_model <- add_criterion(male_fitness_model, criterion = "loo", file = "Fits/male_fitness.model")

loo(male_fitness_model)

## 
## Computed from 8000 by 360 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1561.9 21.4
## p_loo        23.3  1.9
## looic      3123.8 42.8
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     359   99.7%   2275      
##  (0.5, 0.7]   (ok)         0    0.0%   <NA>      
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   1    0.3%   4000      
## See help('pareto-k-diagnostic') for details.

There is one point having a large effect on the posterior. Upon inspection, this data point is not an unreasonable one and there is no cause to remove it from the dataset. It also does not change the causal effect of inheritance treatment on male fitness.

Conduct the posterior predictive check…

pp_check(male_fitness_model, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Derive estimates from posterior

Get predictions from the model and present them in a Table

# predictions for block 1, with UBI GFP

new_data <- tibble(Inheritance_treatment = male_fitness$Inheritance_treatment) %>% 
  distinct(Inheritance_treatment) %>% 
  mutate(Population = 1,
         Block = 1,
         Rearing_vial = 1,
         GFP = "UBI",
         Total_offspring = 1)

predictions_male <- fitted(male_fitness_model, newdata = new_data)

predictions_male <- cbind(new_data, predictions_male)


table1 <-
  predictions_male %>% 
  select(-c(Population, Rearing_vial, GFP, Total_offspring)) %>%
  arrange(Block) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

draws_m <- as_draws_df(male_fitness_model)

# predictions averaged over mediator variables

draws_male <-
  draws_m  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Inheritance_treatment") %>% 
  rename(prop_focal_offspring = value)

draws_diff_m <- draws_m %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Inheritance_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>%
  mutate(`Female - Control` = `Female-limited` - Control,
         `Male - Control` = `Male-limited` - Control,
         `Female - Male` = `Female-limited` - `Male-limited`) %>% 
  select(`Female - Control`, `Male - Control`, `Female - Male`)

# actual Table 1

draws_male %>% 
  group_by(Inheritance_treatment) %>% 
  summarise(`Estimated prop. of offspring sired` = mean(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Inheritance_treatment	Estimated prop. of offspring sired	2.5%	97.5%
Control	0.68	0.574	0.772
Female-limited	0.762	0.662	0.841
Male-limited	0.702	0.594	0.794

# now find the differences between the control and the sex-limited Evolution_treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

prop_table_male <- 
  draws_m %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Inheritance_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Inheritance_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring sired relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Inheritance treatment` = `difference comparison`) #%>% 
  #pander(split.cell = 40, split.table = Inf, round = 3)

Building Figure 2

# female plots

f_1 <-
  draws_female %>% 
  mutate(Inheritance_treatment = fct_relevel(Inheritance_treatment, "Female-limited", "Control", "Male-limited")) %>% 
  ggplot(aes(Inheritance_treatment, prop_focal_offspring)) + 
  stat_halfeye(aes(fill = Inheritance_treatment), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "mean_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5) + # width indicates the uncertainty intervals: we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  coord_flip() +
  xlab("Inheritance treatment") +
  ylab("Female fitness\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_2 <-
  draws_diff %>% 
  gather(key = parameter, value = `Fitness difference`) %>% 
  as_tibble() %>%   

  ggplot(aes(parameter, `Fitness difference`)) + 
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, point_interval = "mean_qi",
               slab_fill = met.brewer("Hiroshige")[1],
               shape = 21, point_size = 3, stroke = 1.5,
               point_fill = "white") + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Treatment contrast") +
  ylab("Female fitness difference\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

# male plots

f_3 <-
  draws_male %>% 
  mutate(Inheritance_treatment = fct_relevel(Inheritance_treatment, "Female-limited", "Control", "Male-limited")) %>% 
  ggplot(aes(Inheritance_treatment, prop_focal_offspring)) + 
   stat_halfeye(aes(fill = Inheritance_treatment), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "mean_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  coord_flip() +
  xlab("Inheritance treatment") +
  ylab("Male fitness\n(prop. offspring produced)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_4 <-
  draws_diff_m %>%
  gather(key = parameter, value = `Fitness difference`) %>% 
  as_tibble() %>% 
  ggplot(aes(parameter, `Fitness difference`)) + 
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, point_interval = "mean_qi",
               slab_fill = met.brewer("Hiroshige")[1],
               shape = 21, point_size = 3, stroke = 1.5,
               point_fill = "white") + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hokusai3", 3)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-0.2, -0.1, 0, 0.1, 0.2)) +
  xlab("Treatment contrast") +
  ylab("Male fitness difference\n(prop. offspring sired)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())


(f_1 + f_2) /(f_3 + f_4) + 
  plot_annotation(tag_levels = 'a')

Figure 2: a shows the estimated distribution of the mean for female fitness for flies carrying autosomes that had previously experienced unconstrained inheritance (control), female-limited inheritance or male-limited inheritance. b shows the difference contrast in female fitness between each of the three inheritance treatments. This difference is on the proportion scale, where a value of 0.1 indicates that females of a given inheritance treatment produce 10 more offspring per every 100 when cohabiting with bw competitor females. c and d depict the same things as a and b, except for male fitness.